clear all;

function retval = SmithWaterman(x,y,sp,sn,g)
	diag = 1;
	up = 2;
	left = 4;
	
	scores = zeros(length(y)+1,length(x)+1);
	directions = zeros(length(y)+1,length(x)+1);
	
	maxScore = 0;
	prevBranch = [];
	prevI = [];
	prevJ = [];
	
	for j = [2:length(y)+1]
		for i = [2:length(x)+1]
			if (x(i-1) == y(j-1))
				scores(j,i) = scores(j-1,i-1) + sp;
			else
				scores(j,i) = scores(j-1,i-1) + sn;
			endif
			directions(j,i) = diag;
			
			score = scores(j-1,i) + g;
			if (score > scores(j,i))
				scores(j,i) = score;
				directions(j,i) = up;
			elseif (score == scores(j,i))
				directions(j,i) += up;
			endif
			
			score = scores(j,i-1) + g;
			if (score > scores(j,i))
				scores(j,i) = score;
				directions(j,i) = left;
			elseif (score == scores(j,i))
				directions(j,i) += left;
			endif
			
			if (0 > scores(j,i))
				scores(j,i) = 0;
				directions(j,i) = 0;
			endif
			
			if (scores(j,i) > maxScore)
				maxScore = scores(j,i);
				prevBranch = [];
				prevI = [];
				prevJ = [];
				prevBranch(end+1) = 1;
				prevI(end+1) = i;
				prevJ(end+1) = j;
			elseif (scores(j,i) == maxScore)
				prevBranch(end+1) = 1;
				prevI(end+1) = i;
				prevJ(end+1) = j;
			endif
		endfor
	endfor
	
	if (maxScore == 0)
		retval = [];
		return;
	endif
	
	i = prevI(end);
	prevI(end) =[];
	j = prevJ(end);
	prevJ(end) =[];
	branch = prevBranch(end);
	prevBranch(end) = [];
	a = '';
	b = '';
	xAlignments = {};
	yAlignments = {};
	while(true)
		if (scores(j,i) == 0)
			branch = 0;		
			xAlignments(end+1) = a;
			yAlignments(end+1) = b;
		endif
		
		switch(branch)
			case diag
				branch = bitshift(branch,1);
				if (bitand(directions(j,i),diag) > 0)
					prevI(end+1) = i;
					prevJ(end+1) = j;
					prevBranch(end+1) = branch;
					i -= 1;
					j -= 1;
					branch = 1;
					a = [x(i),a];
					b = [y(j),b];
				endif
			case up
				branch = bitshift(branch,1);
				if (bitand(directions(j,i),up) > 0)
					prevI(end+1) = i;
					prevJ(end+1) = j;
					prevBranch(end+1) = branch;
					j -= 1;
					branch = 1;
					a = ['-',a];
					b = [y(j),b];
				endif
			case left
				branch = bitshift(branch,1);
				if (bitand(directions(j,i),left) > 0)
					prevI(end+1) = i;
					prevJ(end+1) = j;
					prevBranch(end+1) = branch;
					i -= 1;
					branch = 1;
					a = [x(i),a];
					b = ['-',b];
				endif
			otherwise
				if (length(prevI) == 0)
					break;
				endif
				i = prevI(end);
				prevI(end) =[];
				j = prevJ(end);
				prevJ(end) =[];
				branch = prevBranch(end);
				prevBranch(end) = [];
				a = a(2:end);
				b = b(2:end);
		endswitch
	endwhile
	retval = [xAlignments; yAlignments];
endfunction
